Reproducing the working environment

library(remotes)

# packages that dont build on nix shell
options(repos = "http://cran.rstudio.org")
myPackages = installed.packages()
cranPackages <- c("nnfor", "dtwclust")
to.install <- setdiff(cranPackages, myPackages[, 1])
if (length(to.install) > 0) install.packages(to.install)

if (!("tswgewrapped" %in% installed.packages())) {
    install_github("josephsdavid/tswgewrapped")
}

library(remotes)
options(repos = "http://cran.rstudio.org")

myPackages = installed.packages()
shell <- readLines("shell.nix")
shell <- shell[grepl(shell, pattern = "rPackages")]
shell <- gsub("rPackages.", "", shell)
qhell <- gsub(" ", "", shell)
otherPackages <- c("nnfor", "dtwclust")
cranPackages <- c(otherPackages, shell)

to.install <- setdiff(cranPackages, myPackages[, 1])
if (length(to.install) > 0) install.packages(to.install)

if (!("tswgewrapped" %in% installed.packages())) {
    install_github("josephsdavid/tswgewrapped")
}

Reproducing the models

To reproduce the models from this project, you have two options: you can either load the saved .Rda files, or run all the models (as shown in this report). Note this will likely take several hours, as some of the models take a long time to train, and multiple grid searches over 10,000+ hyperparameters were performed.

Data Definition

Open to the public is data from 5 cities:

• Beijing

• Chengdu

• Guangzhou

• Shanghai

• Shenyang

Each of these datasets contains:

• Air quality measurements at multiple locations (including a US Embassy post)

• Humity

• Air pressure

• Temperature

• Precipitation

• Wind direction

• Wind Speed

Measured at every hour. For detailed info on the data quality and collection methods, please refer to this publication

Pre-processing pipeline

The next step in our analysis workflow is to get the data ready for analysis. First, we will have to do a few tasks

• import data

• Represent data as a time series

• inspect for NAs

• figure out how to deal with NAs

• store in a conveniently manipulatable format

All of this can be seen in R/preprocessing.R, but we will do an in depth explanation here

library(functional)  # to compose the preprocessing pipeline
library(data.table)  # to read csvs
library(rlist)  # for list manipulations
library(pipeR)  # fast, dumb pipes
library(imputeTS)  # to impute NAs
library(pander)  # so i can read the output
library(foreach)  # go fast
library(doParallel)  # go fast

import <- function(path) {
    # first we list the files in our path, in this case 'data/'
    files <- list.files(path)
    # then we search for csvs
    files <- files[grepl(files, pattern = ".csv")]
    # paste the csv names onto the filepath(data/whatever.csv)
    filepaths <- sapply(files, function(x) paste0(path, x))
    # read them in to a list
    out <- lapply(filepaths, fread)
    # Get rid of .csv in each filename
    fnames <- gsub(".csv", "", files)
    # get rid of the confusing numbers
    fnames <- gsub("[[:digit:]]+", "", fnames)
    # set the names of the list to be the clean filenams
    names(out) <- fnames
    out
}

naCount <- function(xs) {
    rapply(xs, function(x) sum(is.na(x)/length(x)), how = "list")
}
# load in the data now
datadir = "data/"
china <- import(datadir)
pander(naCount(china))

tots <- function(v) {
    ts(v, frequency = 365 * 24)
}

totslist <- function(df) {
    # define the column names which we do not want to alter
    badlist <- c("No", "year", "month", "day", "hour", "season", "cbwd")
    # get the column names of the original data frame
    nms <- colnames(df)
    # coerce it into a list
    df <- as.list(df)
    # if the name of the item of the list is in our to not change list, do
    # nothing otherwise, convert to ts
    for (name in nms) {
        if (name %in% badlist) {
            df[[name]] <- df[[name]]
        } else {
            df[[name]] <- tots(df[[name]])
        }
    }
    # return the created list
    df
}

totsall <- function(xs) {
    lapply(xs, totslist)
}

imp_test <- function(v) {
    out <- try(na.interpolation(v, "spline"))
    ifelse(is.ts(out), return(out), return(v))
}

impute <- function(xs) {
    foreach(i = 1:length(xs), .final = function(x) {
        setNames(x, names(xs))
    }) %dopar% imp_test(xs[[i]])
}

impute_list <- function(xs) {
    lapply(xs, impute)
}

to_hash <- function(xs) {
    list2env(xs, envir = NULL, hash = TRUE)
}

preprocess <- Compose(import, totsall, impute_list, to_hash)

Helper functions

First, we define some helper functions, in order to save us time and make our code more readable. First, lets steal some stuff from forecast and tswge:

# pretty plots for our analysis
seasonplot <- forecast::ggseasonplot
subseriesplot <- forecast::ggsubseriesplot
lagplot <- forecast::gglagplot
sampplot <- tswge::plotts.sample.wge
# clean up errant data points
clean <- forecast::tsclean

Next let us define a function which creates a function to resample time series, and then define ones with new sampling rates. Finally, we will write a function which applies each resampling function to a vector

# resample generator
change_samples <- function(n) {
    function(xs) {
        out <- unname(tapply(xs, (seq_along(xs) - 1)%/%n, sum))
        out <- ts(out, frequency = (8760/n))
        out
    }
}

# daily and weekly sampling, monthly is 4 weeks
to_daily <- change_samples(24)
to_weekly <- change_samples(24 * 7)
to_monthly <- change_samples(24 * 7 * 4)
to_season <- change_samples(24 * (365/4))

# pipelining final cleaning and conversion, removing outlier and the couple
# of errant negative values (which do not make sense)
cleandays <- function(xs) {
    xs %>>% clean %>>% abs %>>% to_daily
}

cleanweeks <- function(xs) {
    xs %>>% clean %>>% abs %>>% to_weekly
}
cleanmonths <- function(xs) {
    xs %>>% clean %>>% abs %>>% to_monthly
}
cleanseas <- function(xs) {
    xs %>>% clean %>>% abs %>>% to_season
}

# resample the whole ts
resample <- function(xs) {
    day <- xs %>>% cleandays %>>% window(end = 6)
    week <- xs %>>% cleanweeks %>>% window(end = 6)
    month <- xs %>>% cleanmonths %>>% window(end = 6)
    seas <- xs %>>% cleanseas %>>% window(end = 6)
    list(day = day, week = week, month = month, season = seas)
}

Wge sample plots

First lets load in and resample our data:

library(ggthemes)
library(ggplot2)
library(cowplot)
china <- preprocess(datadir)

#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric

bj <- china$BeijingPM_$PM_US
bj <- resample(bj)

dy <- sampplot(bj$day)

wy <- sampplot(bj$w)

my <- sampplot(bj$m)

sy <- sampplot(bj$s)

Seasonal Plots

Lets check out this seasonal pattern a bit more with a seasonal plot. We will plot them all at once this time and then discuss:

sda <- seasonplot(bj$day) + scale_color_hc() + theme_hc() + ggtitle("Seasonal plot: Daily")
sdw <- bj$week %>>% seasonplot + theme_hc() + scale_color_hc() + ggtitle("Seasonal plot: Weekly")
sdm <- bj$month %>>% seasonplot + theme_hc() + scale_color_hc() + ggtitle("Seasonal plot: Monthly")
sds <- bj$seas %>>% seasonplot + theme_hc() + scale_color_hc() + ggtitle("Seasonal plot: Quarterly")
plot_grid(sda, sdw, sdm, sds, ncol = 2)

Setup

The setup for this bit of analysis is as follows:

china <- preprocess("data/")

#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric
#> Error in na.interpolation(v, "spline") : Input x is not numeric

bj <- china$BeijingPM_$PM_US
bjs <- resample(bj)
bjus <- bjs$day

# split
train <- window(bjus, end = 5)
test <- window(bjus, start = 5)

# clean up
rm(bj, bjus, bjs)

S3 objects and methods

To make analysis of models from different libraries with different structures easier, I wrote some S3 classes and generic methods in order to view and evaluate them. These are shown below

scores <- function(obj) {
    UseMethod("scores")
}

as.wge <- function(x) structure(x, class = "wge")

ASE <- function(predicted, actual) {
    mean((actual - predicted)^2)
}

MAPE <- function(predicted, actual) {
    100 * mean(abs((actual - predicted)/actual))
}

checkConfint <- function(upper, lower, actual) {
    (actual < lower) | (actual > upper)
}

confScore <- function(upper, lower, actual) {
    rawScores <- ifelse(checkConfint(upper, lower, actual), 1, 0)
    sum(rawScores)
}

scores.wge <- function(xs) {
    mape <- MAPE(xs$f, test)
    ase <- ASE(xs$f, test)
    confs <- confScore(xs$ul, xs$ll, test)
    c(MAPE = mape, ASE = ase, Conf.Score = confs)
}

.testPredPlot <- function(xs) {
    p <- ggplot() + theme_hc() + scale_color_hc()
    doplot <- function(df) {
        p <<- p + geom_line(data = df, aes(x = t, y = ppm, color = type))
    }
    out <- lapply(xs, doplot)
    out[[2]]
}


autoplot.wge <- function(obj) {
    testdf <- data.frame(type = "actual", t = seq_along(test), ppm = as.numeric(test))
    preddf <- data.frame(type = "predicted", t = seq_along(obj$f), ppm = as.numeric(obj$f))
    dfl <- list(testdf, preddf)
    .testPredPlot(dfl)
}

Differencing

Before we begin our analysis, we will create a seasonal and trend based difference of our original time series, mostly for the sake of thouroughness. The difference function used is as follows

difference

#> function (type, x, n) 
#> {
#>     szn_trans <- function(x, n) {
#>         artrans.wge(x, phi.tr = c(rep(0, n - 1), 1))
#>     }
#>     arima_trans <- function(x, n) {
#>         f <- artrans.wge(x, phi.tr = 1)
#>         if (n == 1) {
#>             res <- f
#>             return(res)
#>         }
#>         else {
#>             arima_trans(f, n - 1)
#>         }
#>     }
#>     if (is.character(enexpr(type)) == F) {
#>         type <- as.character(enexpr(type))
#>     }
#>     if (type %in% c("arima", "ARIMA", "Arima")) {
#>         return((arima_trans(x, n)))
#>     }
#>     if (type %in% c("ARUMA", "Aruma", "aruma", "Seasonal", "seasonal")) {
#>         szn_trans(x, n)
#>     }
#> }
#> <bytecode: 0x1e740e88>
#> <environment: namespace:tswgewrapped>

trainSea <- difference(seasonal, train, 365)

trainTrend <- difference(arima, train, 1)

Estimation of order

To estimate the order, first all the train sets were put into a list:

trainers <- list(ARMA = train, Seasonal = trainSea, ARIMA = trainTrend)

Then, a wrapper around aic.wge speed up (especially with parallel computation) aic5.wge, was used on the trainers. The speed differences are minor but still much time is saved. The wrappers are:

tswgewrapped:::expand

#> function (v1, v2) 
#> {
#>     c(sapply(v1, function(x) rep(x, length(v2))))
#> }
#> <bytecode: 0x154a9f38>
#> <environment: namespace:tswgewrapped>

tswgewrapped:::rewrite

#> function (v1, v2) 
#> {
#>     rep(v1, length(v2))
#> }
#> <bytecode: 0xd948eb8>
#> <environment: namespace:tswgewrapped>

tswgewrapped:::getpq

#> function (x, p = 8, q = 5, type = "aic", silent = TRUE) 
#> {
#>     if (silent == FALSE) {
#>         cat("Calculating ", type, " for ARMA(", p, ", ", q, " )\n", 
#>             sep = "")
#>     }
#>     res <- try(aic.wge(x, p, q, type))
#>     if (is.list(res)) {
#>         out <- c(res$p, res$q, res$value)
#>     }
#>     else {
#>         out <- c(p, q, 9999)
#>     }
#>     out
#> }
#> <bytecode: 0x13bccdb8>
#> <environment: namespace:tswgewrapped>

aic5

#> function (x, p = 0:8, q = 0:5, type = "aic", silent = TRUE) 
#> {
#>     ip <- expand(p, q)
#>     iq <- rewrite(q, p)
#>     out <- mapply(function(v1, v2) getpq(x, v1, v2, type, silent), 
#>         ip, iq)
#>     out <- (as.data.frame(t(out)))
#>     colnames(out) <- c("p", "q", type)
#>     head(out[order(out[, 3], decreasing = F), ], 5)
#> }
#> <bytecode: 0xb4c00e8>
#> <environment: namespace:tswgewrapped>

tswgewrapped:::aics

#> [[1]]
#> function (x, p = 0:8, q = 0:5, silent) 
#> aic5(x, p, q, type = "aic", silent)
#> <bytecode: 0xa363cc8>
#> <environment: namespace:tswgewrapped>
#> 
#> [[2]]
#> function (x, p = 0:8, q = 0:5, silent) 
#> aic5(x, p, q, type = "bic", silent)
#> <bytecode: 0xa364390>
#> <environment: namespace:tswgewrapped>

aicbic

#> function (vec, p = 0:8, q = 0:5, parallel = FALSE, cl = NULL, 
#>     silent = FALSE) 
#> {
#>     if (parallel == TRUE) {
#>         parLapply(cl, aics, function(f) f(vec, p, q, silent = TRUE))
#>     }
#>     else {
#>         lapply(aics, function(f) f(vec, p, q, silent))
#>     }
#> }
#> <bytecode: 0xd164060>
#> <environment: namespace:tswgewrapped>

The rest of the source code can be seen here. The calculation of all the AICs and BICs is shown below

aics <- lapply(trainers, aicbic, 0:10, 0:8)

pander(aics)

• ARMA:

    *
  
      -------------------------
           p    q    aic
      -------- ---- --- -------
       **90**   9    8   14.62
  
       **53**   5    7   14.62
  
       **68**   7    4   14.62
  
       **11**   1    1   14.62
  
       **99**   10   8   14.62
      -------------------------
  
    *
  
      ------------------------
           p   q    bic
      -------- --- --- -------
       **11**   1   1   14.63
  
       **19**   2   0   14.63
  
       **4**    0   3   14.63
  
       **20**   2   1   14.63
  
       **12**   1   2   14.63
      ------------------------

• Seasonal:

    *
  
      ------------------------
           p    q   aic
      -------- ---- --- ------
       **50**   5    4   15.3
  
       **51**   5    5   15.3
  
       **72**   7    8   15.3
  
       **70**   7    6   15.3
  
       **99**   10   8   15.3
      ------------------------
  
    *
  
      ------------------------
           p   q    bic
      -------- --- --- -------
       **11**   1   1   15.32
  
       **3**    0   2   15.33
  
       **19**   2   0   15.33
  
       **4**    0   3   15.33
  
       **12**   1   2   15.33
      ------------------------

• ARIMA:

    *
  
      ------------------------
           p   q    aic
      -------- --- --- -------
       **58**   6   3   14.63
  
       **38**   4   1   14.63
  
       **49**   5   3   14.63
  
       **44**   4   7   14.63
  
       **20**   2   1   14.63
      ------------------------
  
    *
  
      ------------------------
           p   q    bic
      -------- --- --- -------
       **20**   2   1   14.64
  
       **4**    0   3   14.64
  
       **38**   4   1   14.65
  
       **56**   6   1   14.66
  
       **49**   5   3   14.66
      ------------------------

Estimation of Parameters

To estimate the parameters, we will be using the wrapper estimation function I wrote (just more convenient syntax), as well as the wrapper around ljung.wge I wrote, which is just more convenient to read (and doesnt print 100 residuals)

estimate

#> function (xs, p, q = 0, type = "mle", ...) 
#> {
#>     if (q > 0) {
#>         return(est.arma.wge(xs, p, q, ...))
#>     }
#>     else {
#>         return(est.ar.wge(xs, p, type, ...))
#>     }
#> }
#> <bytecode: 0xd884e60>
#> <environment: namespace:tswgewrapped>

tswgewrapped:::silence

#> function (x) 
#> {
#>     sink("/dev/null")
#>     x
#>     sink()
#> }
#> <bytecode: 0x10d7fb90>
#> <environment: namespace:tswgewrapped>

hush

#> function (f) 
#> {
#>     silence(res <- f)
#>     return(res)
#> }
#> <bytecode: 0x15159038>
#> <environment: namespace:tswgewrapped>

ljung_box

#> function (x, p, q, k_val = c(24, 48)) 
#> {
#>     ljung <- function(k) {
#>         hush(ljung.wge(x = x, p = p, q = q, K = k))
#>     }
#>     sapply(k_val, ljung)
#> }
#> <bytecode: 0x89f8b50>
#> <environment: namespace:tswgewrapped>

All estimations were run with different orders until the best results were found:

estARMA <- estimate(train, 9, 8)

#> 
#> Coefficients of Original polynomial:  
#> 0.1142 0.2027 0.4502 -0.5994 0.0159 0.0679 0.7182 0.0737 -0.1559 
#> 
#> Factor                 Roots                Abs Recip    System Freq 
#> 1-1.3528B+0.9692B^2    0.6979+-0.7381i      0.9845       0.1295
#> 1+1.5505B+0.9575B^2   -0.8097+-0.6235i      0.9785       0.3956
#> 1-0.9753B              1.0253               0.9753       0.0000
#> 1+0.5870B+0.8377B^2   -0.3503+-1.0349i      0.9153       0.3020
#> 1+0.4933B             -2.0270               0.4933       0.5000
#> 1-0.4169B              2.3985               0.4169       0.0000
#>   
#> 

ljung_box(estARMA$res, 9, 8)

#>            [,1]             [,2]            
#> test       "Ljung-Box test" "Ljung-Box test"
#> K          24               48              
#> chi.square 12.9592          31.2323         
#> df         7                31              
#> pval       0.07310876       0.4545512

par(mfrow = c(1, 2))
acf(estARMA$res)
plot(estARMA$res)

estSeas <- estimate(trainSea, 7, 6)

#> 
#> Coefficients of Original polynomial:  
#> -0.1091 -0.1268 0.5919 -0.3199 -0.0725 -0.4779 0.2916 
#> 
#> Factor                 Roots                Abs Recip    System Freq 
#> 1+1.5049B+0.9193B^2   -0.8184+-0.6464i      0.9588       0.3936
#> 1+0.5136B+0.8599B^2   -0.2987+-1.0362i      0.9273       0.2947
#> 1-1.3829B+0.7006B^2    0.9869+-0.6733i      0.8370       0.0953
#> 1-0.5265B              1.8993               0.5265       0.0000
#>   
#> 

ljung_box(estSeas$res, 7, 6)

#>            [,1]             [,2]            
#> test       "Ljung-Box test" "Ljung-Box test"
#> K          24               48              
#> chi.square 18.49173         44.33813        
#> df         11               35              
#> pval       0.07084993       0.1339141

par(mfrow = c(1, 2))
acf(estARMA$res)
plot(estARMA$res)

estTrend <- estimate(trainTrend, 6, 3)

#> 
#> Coefficients of Original polynomial:  
#> -0.0335 -0.6124 0.4784 -0.1861 0.0499 -0.0888 
#> 
#> Factor                 Roots                Abs Recip    System Freq 
#> 1+0.6325B+0.7674B^2   -0.4121+-1.0646i      0.8760       0.3088
#> 1-1.0402B+0.3712B^2    1.4011+-0.8549i      0.6093       0.0872
#> 1+0.4412B+0.3116B^2   -0.7078+-1.6456i      0.5582       0.3147
#>   
#> 

ljung_box(estTrend$res, 6, 3)

#>            [,1]             [,2]            
#> test       "Ljung-Box test" "Ljung-Box test"
#> K          24               48              
#> chi.square 17.52443         43.90103        
#> df         15               39              
#> pval       0.2884916        0.2715692

par(mfrow = c(1, 2))
acf(estARMA$res)
plot(estARMA$res)

Forecasting

We will now proceed with our daily, one year ahead forecasts, as well as the evaluation of them. The fcst function from tswgewrapped is again just a wrappper for fore.whatever.wge:

tswgewrapped::fcst

#> function (type, ...) 
#> {
#>     phrase <- paste0("fore.", enexpr(type), ".wge")
#>     func <- parse_expr(phrase)
#>     eval(expr((!!func)(...)))
#> }
#> <bytecode: 0x11b65758>
#> <environment: namespace:tswgewrapped>

armaCast <- fcst(type = aruma, x = train, phi = estARMA$phi, theta = estARMA$theta, 
    n.ahead = length(test)) %>>% as.wge

seaCast <- fcst(type = aruma, x = train, phi = estSeas$phi, theta = estSeas$theta, 
    s = 365, n.ahead = length(test)) %>>% as.wge

trendCast <- fcst(type = aruma, x = train, phi = estTrend$phi, theta = estTrend$theta, 
    d = 1, n.ahead = length(test)) %>>% as.wge

casts <- list(arma = armaCast, seasonal = seaCast, arima = trendCast)
data.frame(list.rbind(lapply(casts, scores)))

autoplot(armaCast) + ggtitle("ARMA")
autoplot(seaCast) + ggtitle("Seasonal")
autoplot(trendCast) + ggtitle("ARIMA")

Discussion

Although our forecasts were not very good, especially our ARMA forecast will provide us with a nice benchmark to compare our newer models to, at least by ASE. I think the seasonality is an important factor in this dataset, and although our seasonal model did not do so well, it is worth looking more into. We will proceed next with some multiseasonal approaches. As it is daily data, there tend to be two major seasonalities: a weekly pattern, and a seasonal pattern.

S3 methods

We need to construct a few S3 methods for this section. First, the forecast function for Arima objects (those constructed by auto.arima and Arima) returns a forecast of the same length as the length of the external regressors (regardless of the value of h). To deal with this, we write a new method called newFore, which creates a new Arima object that can make forecasts off the length of the test set, and then does a forecast using those xregressors. This returns the same value as running:

predict(`Arima object`, newxreg = `regressors from the test set`, n.ahead = whatever)

but gives us the benefits of the forecast object, such as nice plots, confidence intervals, etc.

newFore <- function(...) {
    UseMethod("newFore")
}
newFore.Arima <- function(obj, newdata, xreg = NULL, h = 1) {
    refit <- Arima(newdata, model = obj, xreg = xreg)
    forecast(refit, h = h, xreg = xreg)
}

Next we want to add new methods for these forecast objects, so we can view our results in a consitent manner:

as.fore <- function(x) structure(x, class = "fore")

autoplot.fore <- function(obj) {
    testdf <- data.frame(type = "actual", t = seq_along(test), ppm = as.numeric(test))
    preddf <- data.frame(type = "predicted", t = seq_along(test), ppm = as.numeric(obj$fitted[1:length(test)]))
    dfl <- list(testdf, preddf)
    .testPredPlot(dfl)
}

scores.fore <- function(obj) {
    mape <- MAPE(obj$fitted[1:length(test)], test)
    ase <- ASE(obj$fitted, test)
    confs <- confScore(obj$upper[1:length(test)], obj$lower[1:length(test)], 
        test)
    c(MAPE = mape, ASE = ase, Conf.Score = confs)
}

Building the model

First we construct a fourier expansion of both the train and test set:

trainExp <- fourier(train, K = c(3, 100))
testExp <- fourier(test, K = c(3, 100))

Next, we plug this in to the auto.arima function, and build our model. Note we are passing the parameter lambda = 0 into the function. Internally, auto.arima does a Box-Cox transformation on the data before running it, to make the model work smoothly. In our case, we do not want values less than zero, because that wouldnt make sense. So, we set lambda, the parameter passed into the box cox transformation to 0, which represents a logarithmic transformation. In turn, this will keep our model, and thus our predictions

mseadyn <- auto.arima(train, seasonal = FALSE, lambda = 0, xreg = trainExp)

Lets check out where the roots lie on the unit circle:

autoplot(mseadyn)

And lets see what order model it identified within the residuals:

data.frame(p = length(mseadyn$model$phi), q = (mseadyn$model$theta), d = length(mseadyn$model$Delta))

Forecasting

First, lets look at what the origina forecast looks like, without doing the newFore function:

forecast(mseadyn, xreg = trainExp) %>>% autoplot

Whew, 5 years ahead. Thats a bit ridiculous. Now lets run the newFore on it, and then check it out again

mseafor <- newFore(mseadyn, newdata = test, xreg = testExp, h = 1)
autoplot(mseafor)

Note the scale on this plot it is now a one year ahead forecast, which is exactly what we want. Unlike the basic ARUMA forecast we did earlier, the confidence limits also lie above zero in all cases, and hopefully it will be a bit better.

Model Evaluation

This model seems to be slightly more appropriate than the simpler models we fit, as it is highly likely this data is multiseasonal, and I much prefer having a lower ARMA order than (9,8). Lets check it out

mseafor %>>% as.fore %>>% autoplot

mseafor %>>% as.fore %>>% scores %>>% t %>>% data.frame

S3 methods

tbats comes from the forecast package, and has similar issues with forecasting as we saw above (its forecasts are the same length as the data it was trained on, h is ignored). So lets write a newFore method for bats objects

newFore.bats <- function(obj, newdata, h = 1) {
    refit <- tbats(model = obj, y = newdata)
    forecast(refit, h)
}

Thankfully, the structure of a tbats forecast is exactly the same as the structure of our dynamic regression forecasts, so for the rest of the methods, we can just coerce our tbats forecasts into fore objects.

Model Building

First, lets note whether train and test are msts or ts objects, just for sanity

sapply(list(train, test), class)

#>      [,1]   [,2]  
#> [1,] "msts" "msts"
#> [2,] "ts"   "ts"

Ok great we are good to go.

cores <- parallel::detectCores()
# > [1] 12
bjbats <- tbats(train, use.parallel = TRUE, num.cores = cores - 1)
autoplot(bjbats)

This is interesting too, so we see the I think trend that the algorithm smoothed out of the time series. Lets go ahead and make the forecast and see how we did

load("analysis/daily/tbats.Rda")

Forecasting

Since we already know what will happen if we just make a forecast with the forecast function, lets instead use our newFore method without comparison

batF <- newFore(obj = bjbats, test, h = 366)
autoplot(batF)

Interesting. What is the autoplot showing us? It is a good question. It is hard to say whether this model is any good or not just from this. Lets try looking at the actual forecasted values with our autoplot method

Model Evaluation

forbats <- as.fore(batF)
autoplot(forbats)

scores(forbats) %>>% data.frame

S3 methods

We will give the usual S3 methods. Because prophet forecasts include the entire dataset, we will have to write methods that take out the training parts. First, we define the coercion method

as.proph <- function(x) structure(x, class = "proph")

Next we define our autoplot and scoring methods:

scores.proph <- function(obj) {
    c(ase = ASE(obj$yhat[-(1:length(train))], test), mape = MAPE(obj$yhat[-(1:length(train))], 
        test))
}

autoplot.proph <- function(obj) {
    testdf <- data.frame(type = "actual", t = seq_along(test), ppm = as.numeric(test))
    preddf <- data.frame(type = "predicted", t = seq_along(test), ppm = as.numeric(obj$yhat[-(1:length(train))]))
    dfl <- list(testdf, preddf)
    .testPredPlot(dfl)
}

Building a model

First, we need to load the data in the proper format

st <- as.Date("2010-1-1")
en <- as.Date("2013-12-31")
trainDates <- seq(st, en, by = "day")
# data MUST BE IN THIS EXACT FORMAT
traindf <- data.frame(ds = trainDates, y = as.numeric(train))

Next, lets make our model. This is pretty easy

library(prophet)
model <- prophet(traindf)

Forecasting

First, we must make a data frame for the future, because prophet has a silly API

future <- make_future_dataframe(model, periods = 366)

Now we can forecast, with the predict function

fore <- predict(model, future)
plot(model, fore)

We can also view how it modeled/forecasted all the components

prophet_plot_components(model, fore) + theme_hc()

#> NULL

Model Evaluation

Lets check out how we did. First, lets take a look at our forecast:

pfore <- as.proph(fore)
autoplot(pfore)

Functions

First, we are now going to have to convert our non time series/numerical data to daily, which means we will have to find a new method of resampling them. For that, we are going to find the most common categorical value in each day, and then take that. We will use the .mode function, defined below:

.mode <- function(x) {
    ux <- unique(x)
    ux[which.max(tabulate(match(x, ux)))]
}

.mode(c(rep("cat", 4), rep("dog", 8), rep("moose", 1)))

#> [1] "dog"

Next we are going to write a function which lumps the categorical variables into groups of each day, takes the .mode of each group, and returns a new vector in that format.

For the lumping, we will use the %/% operator, which returns an integer specifying how many times a number divides another number. For instance:

((1:50) - 1)%/%5

#>  [1] 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6
#> [36] 7 7 7 7 7 8 8 8 8 8 9 9 9 9 9

Then, we will use the glorious tapply function, which does the same thing as the apply family, except it applies the function onto a grouped variable (table apply).

.dailyMode <- function(xs) {
    out <- tapply(xs, (seq_along(xs) - 1)%/%24, .mode)
    out <- unname(out)
    out[1:1826]  # the number of days in 6 years, because we cleaned away extra days
}

Next we are going to write a function that applies this (and our daily time series conversion) onto a list of time series objects and categorical, and spit out a data frame:

dlist <- function(xs) {
    tsbj <- xs %>>% purrr::keep(is.ts)
    tsnobj <- xs %>>% purrr::discard(is.ts)
    out <- append(lapply(tsbj, function(x) window(cleandays(x), end = 6)), lapply(tsnobj, 
        .dailyMode))
    as.data.frame(out)
}

Finally, we are going to write a function that creates a train set and a test set of data frames of time series and categorical variables. Note the use of the <<- superassignment operator. This causes the value to be stored one level out, so if it is a single level function, this will cause the value to be stored in the global environment:

dfsplit <- function(df, n = 5) {
    index <- length(window(df[[1]], end = n))
    tsdf <- df %>% purrr::keep(is.ts)
    nodf <- df %>% purrr::discard(is.ts)
    outs <- lapply(tsdf, function(x) window(x, end = n))
    outn <- nodf[1:index, ]
    train <<- as.data.frame(append(outs, outn))
    tests <- lapply(tsdf, function(x) window(x, start = n))
    test <<- as.data.frame(append(tests, nodf[(index):nrow(nodf), ]))
}

Note we have to do things in this weird way, because otherwise it is extremely difficult to preserve the class ts, which we would like for time series data to be in. Just cbinding or anything like that breaks it immediately.

Analysis of wind direction:

Lets see if the wind direction is of any use to us:

ggplot(train) + geom_point(aes(x = No, y = PM_US.Post, color = cbwd)) + geom_line(aes(x = No, 
    y = PM_US.Post, alpha = 0.1)) + geom_smooth(aes(x = No, y = PM_US.Post), 
    method = "auto") + theme_hc() + scale_color_hc()

Just upon a glimpse this is not that useful. Lets look at a table:

train %>% arrange(cbwd) %>% group_by(cbwd) %>% summarise(mean = mean(PM_US.Post))

It looks like it is an interesting factor to some degree, but at this scale with our mode technique, it is probably strongly distorted. Still it is worth a look. We will note this and examine it when we are performing our models.

Analysis of all numeric variables

plot_vs_response <- function(x) {
    plot(train$PM_US.Post ~ train[[x]], xlab = x)
    lw1 <- loess(train$PM_US.Post ~ train[[x]])
    j <- order(train[[x]])
    lines(train[[x]][j], lw1$fitted[j], col = "red", lwd = 3)
}


numNames <- train %>% purrr::keep(is.numeric) %>% names
numNames <- numNames[-c(1:4, 11:17)]

par(mfrow = c(2, 3))
lapply(numNames, plot_vs_response)

So all of these appear to be fairly interesting i think. However, the summing technique we performed does not make these plots easily interpretable. Mathematically however, they likely still work out, as it is the relationship that matters. When we do linear regression, whether or not we have the total temperature sum or the mean (sum/n) should not make a difference in the significance of the relationships, only the slopes. This means we should be ok. Precipitation does not look useful.

CCF analysis

Lets look at the cross correlation between all the variables next:

ppm <- train$PM_US.Post
ccfplot <- function(x) {
    ccf(ppm, train[[x]], main = x)
}
ccfs <- lapply(numNames, ccfplot)

It appears that low lags are where these variables have the most cross correlation, other than temperature, which seems to have a bigger CCF at a pretty far from zero lag.

Next, to further our EDA, lets perform a cluster analysis of the time series data

More info on clustering time series

source

Moving Forward

Now that we have a basic understanding of how time series clustering works, we can try it out.

Clustering setup

library(dtwclust)

traints <- as.list(keep(train, is.ts))

Now that we have our data and libraries loaded, lets set ourselves up for parallel computation:

# alter for the appropriate number of cores for you my computer has 12, but
# i have a lot of tabs open define the cluster
workers <- makeCluster(10L)
# load the necessary libraries onto the cluster
invisible(clusterEvalQ(workers, library(dtwclust)))
# register the cluster for parallel computation
registerDoParallel(workers)

Partitional Clustering

As we are new to this, lets set ourselves up for a grid search of a bunch of clusters. We will first check out a partitional cluster. For this, we want to use the dtw distance metric, but also check out the DBA centroid. We will also play with their normalization methods, and the time window they use in alignment. we will also try out 2-6 clusters

cfg <- compare_clusterings_configs(types = "partitional", k = 2L:6L, controls = list(partitional = partitional_control(iter.max = 40L)), 
    distances = pdc_configs("distance", partitional = list(dtw_basic = list(window.size = seq(from = 5L, 
        to = 50L, by = 5L), norm = c("L1", "L2")))), centroids = pdc_configs("centroid", 
        share.config = c("p"), dba = list(window.size = seq(from = 5L, to = 50L, 
            by = 5L), norm = c("L1", "L2"))), no.expand = c("window.size", "norm"))

evaluators <- cvi_evaluators("Sil")

comparison <- compare_clusterings(traints, types = "partitional", configs = cfg, 
    seed = 666L, score.clus = evaluators$score, pick.clus = evaluators$pick)

stopCluster(cl)
registerDoSEQ()

comparison$results$partitional %>% arrange(desc(Sil))

clusters <- repeat_clustering(traints, comparison, comparison$pick$config_id)
plot(clusters)

clusts <- clusters@cluster
names(clusts) <- names(traints)
data.frame((clusts))

Hierarchical clustering

Lets try out a bunch of different hierarchical clusters. We should check out different preprocessing methods, different control methods, all the possible centroids, GAK, DTW, SBD, and different window sizes for DTW.

cfg2 <- compare_clusterings_configs(types = "h", k = 2L:6L, controls = list(hierarchical = hierarchical_control(method = "all")), 
    preprocs = pdc_configs("preproc", none = list(), zscore = list()), centroid = pdc_configs("centroid", 
        shape_extraction = list(), default = list(), dba = list()), distances = pdc_configs("distance", 
        hierarchical = list(dtw_basic = list(window.size = seq(from = 5L, to = 50L, 
            by = 5L), norm = c("L1", "L2")), gak = list(), sbd = list())), no.expand = c("window.size", 
        "norm"))

evaluators <- cvi_evaluators("Sil")

comparison2 <- compare_clusterings(traints, types = "hierarchical", configs = cfg2, 
    seed = 666L, score.clus = evaluators$score, pick.clus = evaluators$pick)
comparison2$results$hiera %>% arrange(desc(Sil))

comparison2$results$hierarchical %>% arrange(desc(Sil), distance) %>% ggplot(aes(fill = distance, 
    y = Sil, x = distance)) + geom_boxplot() + facet_wrap(centroid ~ .) + scale_fill_hc() + 
    theme_hc()

cluster2 <- repeat_clustering(traints, comparison2, comparison2$pick$config_id)
plot(cluster2)

plot(cluster2, type = "sc")

Discussion of results

With our clustering, we discovered a few interesting things. Firstly, it appears that by all metrics,air pressure is really different from the rest of the series. Whether that is a good thing or bad thing remains to be seen, but regardless, it is interesting. We will have to pay attention to it as a predictor.

We also need to note how strongly humidity is clustered with our PM2.5 content.This hopefully suggests it will be another strong predictor variable. However, it is important to keep in mind that all our time series are clustered pretty closely, the largest distance is only 0.015. So our clustering may not have had any important results, but at the very least we learned something.

We will now proceed with our multivariate analysis.

Discussion of Appropriateness

In this section we will discuss how each of our weather patterns interplay (or dont) with our air quality.

S3 methods

We need three S3 methods here, an as method, an autoplot method, and a scores method:

as.var <- function(x) structure(x, class = "var")

autoplot.var <- function(obj) {
    testdf <- data.frame(type = "actual", t = seq_along(test$PM_US.Post), ppm = as.numeric(test$PM_US.Post))
    preddf <- data.frame(type = "predicted", t = seq_along(test$PM_US.Post), 
        ppm = (obj$fcst$PM_US.Post[, 1]))
    dfl <- list(testdf, preddf)
    .testPredPlot(dfl)
}

scores.var <- function(obj) {
    mape <- MAPE(obj$fcst$PM[, 1], test$PM_US)
    ase <- ASE(obj$fcst$PM[, 1], test$PM_US)
    c(MAPE = mape, ASE = ase)
}

Order Selection:

First, lets select the possible orders for our VAR forecast

train[1:3] <- NULL
train[9:14] <- NULL
test[1:3] <- NULL
test[9:14] <- NULL
traints <- (purrr::keep(train, is.ts))
traints <- (purrr::keep(train, is.ts))
library(vars)

# simple order
ord <- VARselect(traints, lag.max = 20, type = c("both"))

# Seaosnal order
ords <- VARselect(traints, lag.max = 20, type = "both", season = c(7, 365))

# Another order
ordn <- VARselect(traints, lag.max = 20)

Lets check out what the orders are:

data.frame(rbind(ord$selection, ords$selection, ordn$selection), row.names = c("ord", 
    "ords", "ordn"))

Model Tuning

Next, lets pick the proper parameters with a grid search. First, lets build the grid:

types <- c("const", "trend", "both", "none")
ics <- c("AIC", "HQ", "SC", "FPE")
ps <- c(8, 3, 1)
szn <- c(NULL, 365)
init <- list(ps = ps, types = types, ics = ics, szn = szn)
grid <- expand.grid(init, stringsAsFactors = FALSE)

Then lets set ourselves up for parallel processing:

workers <- makeCluster(11L)
invisible(clusterEvalQ(workers, library(vars)))
registerDoParallel(workers)

Finally, lets run the grid search

gridSearch <- foreach(i = 1:nrow(grid), .combine = rbind) %dopar% {
    model <- VAR(traints, p = grid[i, 1], type = grid[i, 2], ic = grid[i, 3], 
        season = grid[i, 4])
    preds <- predict(model, n.ahead = 366)
    pred <- preds$fcst$PM_US.Post[, 1]
    ASE <- mean((test$PM_US.Post - pred)^2)
    c(ASE = ASE, p = grid[i, 1], type = grid[i, 2], ic = grid[i, 3], season = grid[i, 
        4])
}
stopCluster(workers)
registerDoSEQ()

Lets check out the results:

gridSearch <- gridSearch %>% as.data.frame
gridSearch$ASE <- as.numeric(gridSearch$ASE)
gridSearch[2:4] <- lapply(gridSearch[2:4], factor)
gridSearch %>% arrange(ASE)

Next lets visualize them

library(latticeExtra)
gridV <- (cloud(ASE ~ p + type, gridSearch, panel.3d.cloud = panel.3dbars, scales = list(arrows = F, 
    col = 1), xbase = 0.2, ybase = 0.2, pretty = T, col.facet = level.colors(gridSearch$ASE, 
    at = do.breaks(range(gridSearch$ASE), 20), col.regions = cm.colors), colorkey = list(col = cm.colors, 
    at = do.breaks(range(gridSearch$ASE), 20)), screen = list(z = 40, x = -30)))
gridV

VAR Model: parameter selection and forecasting

We will now go through and select the proper parameters for the model, by calculating the ASE of each VAR prediction, as well as looking at the summary of a linear model (to see significant regressors and help us pick), as well as thinking about our appropriateness assessment:

var1 <- VAR(traints, p = 3, type = "none", ic = "AIC", season = 365)
preds <- predict(var1, n.ahead = 366)
pred <- preds$fcst$PM_US.Post[, 1]
mean((test$PM_US.Post - pred)^2)

#> [1] 3748146

plot(preds)

thelm <- summary(lm(data = traints, PM_US.Post ~ .))

trainsimp <- traints %>% select(PM_US.Post, HUMI, DEWP, TEMP, precipitation, 
    PRES, Iws)
var2 <- VAR(trainsimp, p = 3, type = "none", ic = "AIC", season = 365)
preds <- predict(var2, n.ahead = 366)
pred <- preds$fcst$PM_US.Post[, 1]
mean((test$PM_US.Post - pred)^2)

#> [1] 3747685

trainsimp2 <- traints %>% select(PM_US.Post, HUMI, TEMP, precipitation, PRES, 
    Iws)

var3 <- VAR(trainsimp2, p = 3, type = "none", ic = "AIC", season = 365)
preds <- predict(var3, n.ahead = 366)
pred <- preds$fcst$PM_US.Post[, 1]
mean((test$PM_US.Post - pred)^2)

#> [1] 3746216

trainsimp3 <- traints %>% select(PM_US.Post, HUMI, TEMP, PRES, Iws)

var4 <- VAR(trainsimp3, p = 3, type = "none", ic = "AIC", season = 365)
preds <- predict(var4, n.ahead = 366)
pred <- preds$fcst$PM_US.Post[, 1]
mean((test$PM_US.Post - pred)^2)

#> [1] 3745787

plot(preds)

We have now reached a point where we cannot improve our model any further (trust me i tried). Lets remember these variables for the rest of our multivariate analysis.

Model assessment

vpred <- as.var(preds)
autoplot(vpred)

data.frame(t(scores(vpred)))

Objectively, this model performed worse than the ARMA/mean baseline model. However, visually, it was only off by a bit. It seems to have gotten the period a bit wrong, while it is accurate in scale. It also has a few issues with the shape, but the forecast does not look unreasonable. This may be an interesting model in our bagged results.

What does the nnetar function do?

Basically, the nnetar function takes in your data, properly scales it for the neural network. We will discuss scaling extensively in the next section, do not worry. After scaling, if needed, it calculates the AR order of the time series, as well as the lags. Then, it takes your scaled time series, the things it calculated and adjustements it made, into the nnet function from the nnet library. Then, you have a model and you make predictions. It is very user friendly, and very very fast.